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The intergalactic medium was reionized before redshift z ~ 6, most likely by starlight which 
escaped from early galaxies. The very first stars formed when hydrogen molecules (H2) cooled gas 
I/"") ■ inside the smallest galaxies, minihalos of mass between 10 5 and 10 8 M©. Although the very first 

\Q I stars began forming inside these minihalos before redshift z ~ 40, their contribution has, to date, 

been ignored in large-scale simulations of this cosmic reionization. Here we report results from the 
first reionization simulations to include these first stars and the radiative feedback that limited 
their formation, in a volume large enough to follow the crucial spatial variations that influenced 
^ ' the process and its observability. We show that, while minihalo stars stopped far short of fully 

'k^J . ionizing the universe, reionization began much earlier with minihalo sources than without, and was 

greatly extended, which boosts the intergalactic electron-scattering optical depth and the large- 
angle polarization fluctuations of the cosmic microwave background significantly. Although within 
current WMAP uncertainties, this boost should be readily detectable by Planck. If reionization 
ended as late as z ov < 7, as suggested by other observations, Planck will thereby see the signature 
of the first stars at high redshift, currently undetectable by other probes. 

Subject headings: cosmology: theory — galaxies: high-redshift — radiative transfer 
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1. Introduction 

The theory of reionization has not yet advanced 
to the point of establishing unambiguously its tim- 
ing and the relative contributions to it from galax- 
ies of different masses. In a Cold Dark Matter 
("CDM") universe, these early galactic sources 
can be categorized by their host halo mass into 
minihalos ("MHs") and "atomic-cooling" halos 
("ACHs"). MHs have masses M ~ 10 5 - 10 s M 
and virial temperatures T v i r ~ 10 4 K, and thus 
molecular hydrogen (H 2 ) was necessary to cool the 
gas below this virial temperature to begin star for- 
mation. ACHs have M > 10 s M Q and T vir > 10 4 
K, for which H-atom radiative line cooling alone 
was sufficient to support star formation. The 
ACHs can be split further into low-mass atomic- 
cooling halos ( "LMACHs" ; M ~ 10 8 - 10 9 M ), 
for which the gas pressure of the photoionization- 
heated intergalactic medium ("IGM") in an ion- 
ized patch prevented the halo from capturing the 
gas it needed to form stars, and high-mass atomic- 
cooling halos ( "HMACHs" ; M > 10 9 M Q ), for 
which gravity was strong enough to overcome this 
"Jeans-mass filter" and form stars even in the ion- 
ized patches. 

Once starlight escaped from galactic halos into 
the IGM to reionize it, the ionized patches ("H 
II regions") of the IGM became places in which 
star formation was suppressed in both MHs and 
LMACHs. At the same time, UV starlight at en- 
ergies in the range 11.2 - 13.6 eV also escaped from 
the halos, capable of destroying the H 2 molecules 
inside MHs through Lyman- Werner band ( "LW" ) 
dissociation, even in the neutral zones of the IGM. 
This dissociation eventually prevented further star 
formation in some of the MHs where the back- 
ground intensity was high enough. Early esti- 
mates, in fact, suggested that this would have 
made the MH contr ibution to reionization small 
200fj ). and, until now, large-scale 



(Haiman et al 



simulations of reionization have neglected them al- 
together. 

In this letter we report the first radiative trans- 
fer (RT) simulations of reionization to include all 
three of the mass categories of reionization source 
halos, along with their radiative suppression, in 
a simulation volume large enough to capture both 
the global mean ionization history and the observ- 
able consequences of its evolving patchiness in a 



statistically meaningful waj0. We overcame the 
limitation of previous large- volume simulations by 
applying a newly developed sub-grid treatment to 
include MH sources (section [5]), and calculating 
the transfer of LW-band radiation self-consistcntly 
with the source population using the scheme of 
~ (|2009h . 



Ahn et al 



2. Methods 

We performed a cosmological N-body simula- 
tion of structure formation with 3072 3 particles in 
a 114/ft.Mpc simulatio n box, using the WM AP5 
background cosmology (jDunklev et al.ll2008 ). For 
this we used the code CubeP 3 M (|Merz et al.ll2005t 
Iliev et all [20081 : J. Harnois-Deraps et al. 2012, 
in preparation), in which the gravity is com- 
puted by a P 3 M (particle-particle-particle-mesh) 
scheme. The simulation was started at redshift 
z = 300 and run to z = 6. N-body data 
were recorded at 86 equally-spaced times (every 
11.53 Myrs) from z = 50 to z — 6. Each data 
time-slice was then used to create matter density 
fields by smoothing the particle data adaptively 
onto a uniform mesh - or an "RT grid" - of 256 3 
cells. All cosmological halos with M > 10 8 M Q 
(corresponding to 20 particles or more), and thus 
both LMACHs and HMACHs, were identified on- 
thc-fly using a spherical overdensity halo finder 
with overdensity of A = 178 with respect to the 
mean. 

Because MHs are too small to be resolved in our 
simulation box, our RT grid was populated with 
MHs through a newly-developed sub-grid model, 
as follows. We started with a separate, high- 
resolution N-body simulation of structure forma- 
tion in a box with (6.3/h Mpc) 3 volume and 
1728 3 particles, which resolved all MHs with M > 
10 5 M Q with 20 particles or more. We then par- 
titioned the box into a uniform grid of 14 3 cells, 
such that each cell is the same size as one of the 
RT grid cells in our main, 114/ft, Mpc simulation 
box, and calculated cell density and the total num- 
ber of MHs per cell. A strong and tight correla- 
tion between the number of MHs located in a cell 
and its density is observed (Fig. [I]). The best fit 
to this correlation at each redshift was then used 
as the total number of MHs in each grid cell of 



1 Fi rst results of thes e simulations were briefly summarized 
in I Ahn et all teOloT) . 
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z=10.1; V ceD =(0.64 Mpc) 3 z=20.1; V oell =(0.64 Mpc)= _ _ z = 30.0; V cell =(0.64 Mpc)= 




log 10 (l+<5) log 10 (l + 5) log 10 (l + 5) 



Fig. 1. — Correlation between the total number of MHs per RT cell (Ns-.s) and the cell density in units 
of the mean density (1 + S), based on a 6.3//iMpc box N-body simulation which resolves all halos with 
M > 10 5 M Q . Plots are for correlations at three different redshifts, z =30, 20.1 and 10.1 from left to right. 
The volume of the RT cell is (0.64Mpc) 3 . 



114//i Mpc box. 

Based on these structure formation results for 
the IGM density field and the source dark matter 
halos, we then calculated the radiative transfer of 
H-ionizing and H2-dissociating photons (see Ta- 
ble [T] for the RT simulation parameters). The 
stars inside ACHs are assumed to produce g 1 ion- 
izing photons per baryon every 10 Myrs, where 
9i = /y/CWiOMyr), and where f 1 = f e f*Ni, f c 
is the escape fraction of ionizing photons, /* is the 
star formation efficiency, and Ni is the number of 
ionizing photons per stellar baryon produced over 
the star's lifetime t* - we use t* = 11.53 Myr for 
both HMACHs and LMACHs, and t* = 1.92 Myr 
for MHs. We assign one Pop III star per MH, 
motivated by numerical simulations of first star 
formation inside MHs, which find that typically 
one Pop III star with a mass between 100 and 
1000 Mq forms per MH in t he absence of strong 
soft UV radiative feedback ( Bromm et al.l 12002 ; 
Abel et all 120021 lYoshida et alj l2006h . At each 
cell, only those MHs which are newly collapsed 
every 1.92 Myrs are assumed to host Pop III stars 
to roughly approximate the disruptive radiative 
and mechanical feedback by the first star and its 
by-products (such as a supernova) on the halo gas 
(for this, we create "morphed" density fields every 



1.92 Myrs by linearly interpolating the N-body 
density fields in time which are separated by a 
time interval of 11.53 Myrs, and finite-difference 
corresponding minihalo populations on each cell). 
Star formation in MHs is further suppressed when 
the local LW background calculated at each time 
step in 3D using the scheme bv lAhn et al. (2009), 
but now considering both ACHs and MHs and also 
improved in speed using the fast Fourier trans- 
form (FFT) scheme - reaches a certain thresh- 
old JLw.th- At present the precise value of this 
threshold is not well determined, but the typi- 
cal values found by high-resolution simulations of 
MH star formation are Jlw. th = [0-01 — 0. 11 x 



io- 



erg s cm sr (|Ma chace k et al 



Yoshida et al. 2003; O'Shea fc Norman 



2001 
20081) 



We adopted a constant value chosen from this 
range for each simulation. Even though stars 
may still form by H 2 cooling, when Jlw > 
Jlw, th, in MHs in mass range M ~ 2 x 10 6 — 
10 8 M Q and M ~ 10 7 - 10 s M© at J LW ,th ^ 
0.01 x 10~ 21 ergs -1 cm~ 2 sr -1 and JLw.th — 

-21 , 



0.1 x 10 21 ergs cm sr respectively (e.g. 
O'Shea fc Norman] 120081 ). we neglect their con- 



tribution because they constitute only a small 
fraction of the whole MH population. 

The simulations with ACHs only (and without 
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LW radiative transfer) are described in Iliev et al 
d2012h . Simulation parameters are given in Ta- 
bled] 

3. Role of the first stars during cosmic 
reionization 

We demonstrate the effects of the first stars by 
direct comparison of the results from two sim- 
ulations, a fiducial case which includes all ion- 
izing sources down to the first stars hosted by 
MHs (Case L2M1J1) and a corresponding refer- 
ence case which includes the larger, atomically- 
cooling halos with exactly the same properties, but 
no MH sources (Case L2, previously presented in 



Iliev et alJl2012h . Our results show that the early 



reionization history is completely dominated by 
the first stars, while the late (redshift z < 10) 
history is driven by the stars inside HMACHs 
(Figs[2j\ and [2)3, top panel). The very first stars 
start to form inside MHs at redshift z ~ 40, and 
dominate the reionization process until z ~ 10 but 
through self-regulation which slows their contri- 
bution to reionizatiord (Fig. [2)3). Although the 
abundance of ACHs rises exponentially, they re- 
main relatively rare, and thus sub-dominant, un- 
til z ~ 10. After redshift z ~ 8, though, the 
two reionization histories become largely indistin- 
guishable, because the same HMACHs then dom- 
inate reionization and push Jlw above Jlw, th (at 
z ~ 12) halting MH star formation altogether, 
long before the MHs can complete reionization on 
their own. 

Nonetheless, the MH sources (the first stars) 
can have quite a dramatic effect on the electron- 
scattering optical depth r os . While intergalactic 
H II regions fully overlap (at redshift z ov , here 
defined as when the mass-weighted mean ionized 
fraction in the IGM, x m , first surpassed 99%) at 
almost identical redshifts z ov ~ 6.8 with (L2M1J1) 
or without (L2) MHs, the early rise of x m with 
MH sources boosts the optical depth by as much as 
47% relative to that without MH sources: r os = 
0.0861 for L2M1J1, while t cs = 0.0603 for L2. 
This satisfies the current observational constraints 



2 The oscillation of Jlw around the plateau at 
^Lw/^LW,th ~ 1 observed in Figs. [2)3 and \3\ is a 
numerical artifact which occurs because LW suppression 
locks the MH star formation rate onto the level that keeps 
•^LW = ^LW th, an d the simulation time step (1.92 Myrs) 
is comparable to the MH formation time scale. 



on reionization: (1) reionization ended no earlier 
than redshift z = 7 ( Fan fc et al. 20061 [Ota et al 



(Fe 
1. 



2010: iMortlock et al l 
IPentericci et al.ll201llh a nd (2) 



20111 iBolton et all boilF 



.088 ± 0.015 

at 68 % confidence level ([Larson et a,l.ll201lh . Pre- 
dicted values of r es and z ov are model-dependent, 
and thus we tested the robustness of our conclu- 
sions by varying the physical parameters of MHs 
and ACHs (Fig.©. 

The first stars, born inside MHs, imprint a 
distinctive pattern on the global reionization his- 
tory. For example, in case L2M1J1, when the 
LW-platcau ended, reionization briefly stalled, 
since MHs no longer formed the stars which re- 
plenished the ionizing background and only ACH 
sources remained, thereafter; the ACH contri- 
bution took a bit more time to climb enough 
to move reionization forward again. This ex- 
plains the brief "x m -plateau" from z ~ 12 to 
z ~ 10 in Fig. for case L2M1J1, while in 
case L2, x m grows continuously without showing 
such a plateau (Fig. [2j3) . This feature is generic 
(see Figs [3] and [4)3 for different sets of param- 
eters we explored). Reionization histories with- 
out MH sources 
RT si mulat ions 
20071: 



2006 



, modelled either by large-scale 
(ICiardi et all 120031: llliev et al. 
Trac fc Cenl |2007t iMcQuinn et al 



2007: Zahn et al. 2007; Iliev et al. 2012) 



or semi- 



analy tical calculations (jHaiman et al . 2000; Zahn et al 
2007n . are all similar in that respect. 

Reionization histories with MH sources calcu- 
lated here, however, find an ionization plateau 
phase. Previous studies that considered MH 
stars and their impact were not able to settle 
the issue of their global effect on reionization. 
This is either because they simulated volumes 
much too sm all to represent a fair portion of 
the Universe ( Ricotti et ahl 120021: ISokasian et al 
2004; lYoshida et all 120061 : iRicotti et al.l 1200 



Johnson et al.l 20121) . or else treated reionization 
by a semi-analytical, 1-zone, homogeneous ap- 



(Haiman et al. 20001 Furlanetto & Loeb! 2005) 


or without (IShaoiro et al.l 


Il994t IWvithe fc Loeb 



20031 : lHaiman fc BryaiJ200alWyithe fc Cenll2007h . 
which cannot capture its innate spatially inhomo- 
geneous nature, or made a semi-analytical ap- 
proximation that accounted statistically for spa- 
tial inhomogeneity but without LW suppression 



(jKramer et al.l 2006) 
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We find that the global reionization history at 
z < 20 depends on JLW.th more strongly than 
Af*, in- This is due to the very nature of self- 
regulation of the first star formation. The larger 
the Jlw, th is, the weaker the suppression of star 
formation becomes, thereby temporarily hastening 
the progress of re i onization. I f M+ t i t is smaller 



( Turk et all 120091: IStacv et al.l 120101 Weif et al. 
2011a ; Stacv et al. 2012t ) than those simulated 



here, those stars produce less ionizing and LW ra- 
diation and the resulting suppression is weaker, 
which partly compensates for the lower emission 
per star. Similar type of compensation would oc- 
cur also when the number of MHs with a potential 
to form the first stars are smaller than our esti- 
mate, due to the relativ e offset of baryonic gas 
from some of the MHs (jTseliakhovich fe Hirata 
20101 : iGreif et al.ll2011b[ ). 



4. Probing the first stars with Planck 

We thus conclude that the first stars hosted 
by MHs likely made an important contribution to 
reionization. But how can we probe them observa- 
tionally? While significant, the effects of the first 
stars are largely confined to the early stages of 
reionization, at redshifts z > 10, which puts them 
beyond the reach of most current instruments. 
Recent observations by the South Pole Telescope 
(SPT) have been used to place an upper limit on 
the kinetic Sunyaev-Zel'd gvich (kSZ) effect from 
the epoch of reionization (jReichardt et al.l 1201 if) . 
While it has been sugge sted that this rest ricts the 
duration of reionization (jZahn et al.ll201 lh , we will 
show elsewhere (H. Park et al. 2012, in prep) that 
the kSZ signal from our fiducial case, L2M1J1, is 
well below the observed upper bound. The com- 
bined effect of the first stars will be reflected in 
the cosmic microwave background (CMB) polar- 
ization anisotropics at large scales. The current 
best constraints on r cs by the WMAP satellite 
r es are still relatively weak, and thus models with 
low-Tes values like L2 are still acceptable at the 
2cr (95%) confidence level (Fig.[2j3, middle panel). 
However, through the far more precise measure- 
ment of the CMB polarization by the Planck mis- 
sion we should be able to discern the influence of 
the first stars on reionization. 

As we discussed above, current observational 
constraints suggest that reionization was not com- 



plete before z ~ 7. Imposing this condition as a 
prior on the allowed reionization histories x m (z), 
we predict that the Planck mission will clearly de- 
tect the era of first stars (Fig. 2]). In Fig. @] we show 
a statistical measure of the Planck sensitivity to 
detecting the signature of the first stars through 
the pr incipa l component analy s is bv Hu fe Holder 
(|2003l ) and iMortonson fc Hul (120081 w ho modT 
fied COSMOMC (|Lewis fc Bridle! I2002T ) to allow 
generic reionization models Reionization prin- 
cipal components {S^^z)} are eigen-vectors of the 
relevant Fisher matrix (evaluated with an artibrar- 
ily chosen fiducial history x e , &d(z)), 



^ I 1\ d 2 Cf E 



1=2 



2 J dx e (zi)dx e (zj) 



(1) 



which can be used to describe any generic ion- 
ization history x e (z) with just a small number of 
modes, such that 

iVmax 

x e (z) =x e< f Ld (z) + y^^m^S^jz). (2) 

The mode amplitude m M , for a given history x e (z), 
becomes 

J^dzS^z) [x e (z) -x et&d (z)] 
m p = . (3) 

^max ^min 

Based on the Planck data after its full 2 years of 
planned operation, the narrow posterior distribu- 
tion of allowed r es values will allow us to distin- 
guish reionization models like L2 and L2M1J1 un- 
ambiguously, and thereby strongly constrain the 



3 While we use the scheme by Mortonson & Hu (2008), we 
implement the following ingredients to optimize the anal- 
ysis for our purpose. First, to apply the late-reionization 
prior, z ov < 7, we created 7 sets of principal components 
based on x e ^(z) = 4 ^~ z s (see Eq. [2] for the defini- 
tion of x e fidM), which make x e (z) behave well around 
z ~ z ov . We then use this late-reionization prior to re- 
ject any sample reionization history with z ov > 7 when 
forming the Monte-Carlo Markov-chain of varying reion- 
ization models. Second, we improve the physicality con- 
dition, or < x e (z) < 1 at any z, which was somewhat 
poorly applied in Mortonson & Hu (2008). Whenever a 
set of m M parameters (Eq. |3J are sampled, we calculate the 
corresponding x e (z), and when either min(x e (z)) > 0.04 or 
max(a; e (z)) < 1.04 is violated, we reject that sample. This 
small, 4% non-physicality in x e is still necessary because of 
the oscillatory nature of x e (z) caused by the limited num- 
ber of principal components, but has only modest effects 
on the CMB E-mode polarization power spectrum C ; EE . 
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available reionization models. A high measured 
measured value of r cs > 0.085 will be a clear (if 
indirect) signature of the first stars. 

Finally, we note that the presence of MH 
sources introduces the x m -plateau noted above, 
which in turn imprints characteristic features in 
C EE . Hence, Planck might be able to distinguish 
(albeit at lower statistical significance, of > 2a 
or > 95%) reionization models with and without 
first stars even if they have very similar values 
of r es and z ov (Fig. HJ3). Full reionization simu- 
lations like ours find it hard to satisfy both of 
these observational constraints without including 
a significant contribution from t he first stars, but 
some semi-analytical models (|Haiman fe Brv; 
20061 lHaardt fc Madaul 120121 ; #0.348.67.8 and 



#2.609C_165.2 in Fig. |4f3 which are of unnatu- 
rally large gaps in relative efficiencies of LMACH 
and HMACH) do find such scenarios. However, 
all such models lack the plateau feature in x m (z) , 
regardless of the details of the assumed physics, 
and reside in a narrow window of m^-parameter 
space adjacent to that occupied by our no-MH 
cases, as demonstrated in Fig. [4)3. 

In summary, Planck is capable of distinguish- 
ing with high confidence between definitive classes 
of reionization scenarios allowed by the current 
constraints, and thereby significantly restricting 
the available parameter space. Planck will either 
probe the signature of the first stars, or show that 
the first stars had a negligible impact on reion- 
izaion. Once these first results confirm the role of 
the first stars, simulations of the type presented 
here can be used to study other observable quan- 
tities and thus deepen our understanding of the 
early universe. 
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Fig. 2. — (A) Maps of evolving hydrogen-ionized fractions at different rcdshifts (rows), for our fiducial 
model with MH sources included, L2M1 Jf (1st column), vs. the corresponding reference model with only 
atomically-cooling halos, case L2 (2nd column). The slices are 0.45//i Mpc-thick. Color represents linearly- 
scaled ionized fraction from (blue) to 1 (red). (B) (top) Globally-averaged history of the mass- weighted 
ionized fraction for models L2M1J1 (black, solid) and L2 (blue, dashed), (middle) r es integrated from z — 
to redshift z for L2 and L2M1J1. (bottom) Evolution of the mean Jlw hi units of JLw,th for Case L2M1J1. 
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Fig. 3. — Model dependency of the history of cosmic reionization. In addition to cases LI and L2, which 
do not account for MHs, we show predictions for MH-included cases by parametrizing the star formation 
inside MHs through M* t ni (mass of the Pop III star) and Jlw, th (threshold Lyman- Werner intensity). (Left) 
Late-overlap models. (Right) Early-overlap models. 
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Fig. 4.— (A) Detecting the first stars, [left]: Forecasts of Cf E of cases L2M1J1 (with MHs) and L2 (no 
MHs) for Planck. The error bars are es timated Planck 2-year la sensitivity including cosmic variance (top 
panel; iThe Planck Collaboration! l2006f) . [right]: Model-selection power of Planck. Contours represent la 
(68%), 2a (95%) and 3a (99.7%) confidence levels from inside out, on marginalized posterior distributions 
of selected parameters (m M 's and T es ) using mock data based upon Model L2M1J1 (black square). Case L2 
(blue triangle) can be ruled out only from the measurement of T es by Planck. The prior condition of z ov < 7 
is applied, which rules out early reionization (z ov > 8) models. (B) Breaking the degeneracy in z ov and 
r cs . [left]: Ionization histories of various models, but with identical z v(— 6.8) and r cs (~ 0.085). The model 
with MH sources (case L2M1J1, black line) stands out from almost identical, no-MH models. g0.348_67.8 
(no clumping; blue, dotted) and q2.609C 165.2 (with z-dep endent clumping; cyan, dashed) are semi-analytic 
models obtaine d from equation (Al) o f flliev et ah! ( 20071 ) with n — 0.1, and Haardt & Madau (red, dot- 
dashed) is from Haardt &: Madaul f|2012h . [right]: Hypothesis-testing power of Planck on MH-included (black 
square) vs. no-MH models (triangles). Contours have the same meaning as those in (A). No-MH models are 
clustered and well separated from case L2M1J1 at > 2a confidence level. 
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